Quantum Transport in Chemically-modified Two-Dimensional Graphene: From Minimal 

Conductivity to Anderson Localization 
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An efficient computational methodology is used to explore charge transport properties in chemically-modified 
(and randomly disordered) graphene-based materials. The Hamiltonians of various complex forms of graphene 
are constructed using tight-binding models enriched by first-principles calculations. These atomistic models are 
further implemented into a real-space order-N Kubo-Greenwood approach, giving access to the main transport 
length scales (mean free paths, localization lengths) as a function of defect density and charge carrier energy. 
An extensive investigation is performed for epoxide impurities with specific discussions on both the existence 
of a minimum semi-classical conductivity and a crossover between weak to strong localization regime. The 2D 
generalization of the Thouless relationship linking transport length scales is here illustrated based on a realistic 
disorder model. 

PACS numbers: 73.63.-b, 72.15.Lh, 73.63.Fg, 63.22.-m 
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I. INTRODUCTION 

Ever since graphene was experimentally synthesized in 
2004 1 , interest in its promising conduction properties has 
increased considerably^^!. Owing to its two-dimensionality 
and reported large charge mobility, monolayer graphene has 
been initially envisioned as a genuine candidate to replace 
silicon in nano-electronics 4 6 . But despite its realistic poten- 
tial in high-frequency device applications^, the absence of 
a substantial band-gap hinders its use for replacing silicon 
MOSFET devices in logic applications^. 

Various solutions have already been proposed to over- 
come this hurdle such as opening a wide band-gap, us- 
ing quantum confinement in ribbon d 10 " 12 ! or using chemi- 
cal oxidation or hydrogenation to break the symmetry of 
the graphene lattice 13 22 . Both of these met hods h ave how- 
ever been demonstrated to be far too invasive^ 23 ! 24 ! generat- 
ing a large quantity of defects and damaging the otherwise 
Dirac-like properties of electronic excitations. Other more 
seducing proposals include the use of a laser field in the 
mid-infrared range which can induce tunable band gapJSl 
electric -field assisted gap opening in bilayers 26 or chemical 
doping which in certain conditions allow to engineer con- 
trolled mobility gaps as large as 1 e V 27 * 28 ! 

In all cases, the precise understanding of the impact of 
disorder on electronic and (charge, spin and phonon) trans- 
port properties of graphene appears of paramount impor- 
tance. Disorder in graphene exhibits many different fla- 
vors from structural defects to adsorbed impurities, recon- 
structed edges or long range Coulomb scatterers trapped in 
the graphene substrate (oxide layer). To date, the detailed re- 



lationship between microscopic complexity of disorder fea- 
tures and the onset of graphene unique transport properties 
remains elusive. This is particularly debated in relation with 
the so-called Klein tunneling mechanism 2 ^ and the weak 
anti-localization phenomenon which are both manifestations 
of pseudospin effects 31 * 34 . 

Disorder first comes as a source of elastic scattering which 
limits the mean free path in a way which strongly depends 
on the disorder potential characteristics. The energy depen- 
dence of the mean free path and associated semi-classical 
transport quantities such as the Drude conductivity and the 
charge mobility can be indeed connected to the short or long 
range nature of the scattering potential 35 . Beyond the occur- 
rence of a diffusive regime, quantum interferences contribute 
significantly to the transport features at sufficiently low tem- 
peratures. In a dditio n to the conventional weak localiza- 
tion phenomenorP^I, crossovers from weak localization to 
weak anti-locali zation h ave been predicted and experimen- 
tally observed 3 ^ 32 ! 34 ! 38 " 44 ! Pseudospin-related quantum in- 
terferences are however maintained provided disorder does 
not break all underlying symmetries. This is not the case in 
presence of chemical defects which damage the sp 2 lattice 
symmetry. Such stronger disturbances of graphene structure 
maximize localization effects, eventually turning the mate- 
rial to a two-dimensional insulator (Anderson localization). 
If Anderson localization has been highly debated and contro- 
versial for long-range disorder 45 52 , its relevance for strongly 
damaged graphene is now well documented both theoreti- 
cally and experimentally^ 5 ^. A recent theoretical study has 
however related the existence of a robust metallic state in 
presence of local magnetic ordering for partly hydrogenated 
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graphene^2l5°l pinpointing possible subtleties between corre- 
lated impurity distribution and transport features. 

The main objective of this paper is to illustrate how an 
insulating regime can be tuned by intrusive functionaliza- 
tion of a graphene sheet caused by oxygen atoms bound in 
the epoxy position. Epoxide defects are for instance incor- 
porated on graphene after ozone treatmenpSI. These epoxy 
impurities have a drastically different impact on resonant en- 
ergy peaks in the vicinity of the Dirac point when compared 
to single impurities 6 ^. 

To address this objective, the oxygen in epoxy position 
is studied by means of accurate ab initio techniques. This 
model allows us, on the one hand, to prove that the oxygen 
epoxy bonding lies in between a pure sp 3 -like covalent bond 
and an ionic bond, and on the other hand, to supply a suitable 
tight-binding (TB) model for further studies in very large 
scale systems. Using this TB model, the Kubo-Greenwood 
formalism is implemented in real space to obtain meaningful 
transport length scales and conduction properties. Several 
quantities such as the mean free path, the semi -classical con- 
ductivity and the localization length are analyzed in depth. 
The ongoing debate concerning long and short-range scat- 
tering behavior is also briefly commented in light of our re- 
sults. The crossover to the strongly localized regime is then 
investigated. Finally, conventional scaling laws are tested on 
our model in this localization regime. 

II. EPOXY DEFECTS 

V.V. Cheianov et al. llcJTl demonstrated the tendency of 
epoxy-bound adatoms to form spatially correlated states. 
The interaction between epoxy groups is mediated by 
the conduction electrons, similar to the Ruderman-Kittel- 
Kasuya-Yosida (RKKY) interaction which correlates mag- 
netic impurities 62 . These ordered states only exist for low 
impurity densities and disappear at a critical temperature T c . 
For high concentrations of epoxy groups due to oxidatiorP^l 
of graphene, J.L. Li et al. Il64l argue that two epoxy groups 
attaching on the opposite ends of a carbon hexagon create 
more open rings inducing cracks along neighboring rings. 
Similarly, more recently, S. Fabris et al. [65 1 present a mech- 
anism giving rise to more complex crack propagation. How- 
ever, crack propagation has, up to our knowledge, only been 
reported under strong reduction and oxidation treatments. 
Furthermore, an aqueous environment seems to be manda- 
tory. H.J. Xiang et al. [66] also report on these unzipped 
chains caused by epoxy groups thus inducing lower energy 
conformations. Their study is however limited to concentra- 
tions above 25 % of epoxy density. 

The model investigated in the present paper takes ad- 
vantage of this literature while limiting its complexity to 
avoid any loss of generality for the simulations of differ- 
ent moderate concentrations (epoxy density ranging from 
0.01 % to 5 %). Consequently, epoxy groups are assumed 
to be randomly distributed over the graphene sheet and 




FIG. 1: (color online) An oxygen atom in epoxy position on a 5 x 5 
cell (a). 3D (b) and 2D (c) charge density difference as defined in 
the text. Charge accumulation in red/green and charge depletion 
in blue. Isovalues of 0.006 |e|/A 3 and -0.006 je|/A 3 for the 3D 
charge density difference. 



the model prohibits the destructive presence of two oxygen 
atoms on the same hexagon. This simplified model could 
well describe the functionalization of graphene due to ozone 
treatment an d comparisorPl with experimental results^ 
backs this up. 



m. NUMERICAL TECHNIQUES 

The first part of this section presents ab initio calcula- 
tions performed to predict the structural properties of epoxy 
bound oxygen. Likewise, it handles how the TB parameters 
were extracted from these calculations. The second part of 
the section sets out the foundations of the Kubo-Greenwood 
formalism developed in a TB framework. 



A. From ab initio to tight-binding models 

The Density Functional Theo ry (D FT) calculations are 
conducted using the SIESTA cod d 68 " 70 ! within the local den- 
sity approximation (LDA) on the exchange-correlation func- 
tional in the Ceperley- Alder 71 form parametrized by Perdew 
and Zunger 72 . Core electrons are included using Troullier- 
Martins 73 pseudopotentials. Double C, plus polarization or- 
bitals are used to define the basis set. 

Fig.[TJa) illustrates the oxygen atom in epoxy position af- 
ter ab initio geometry optimization. In order to investigate 
the chemical bonding of oxygen on graphene, the charge 
density difference is calculated. Plotting such charge den- 
sity difference allows for a visual understanding on how the 
electronic clouds associated with the orbitals are altered by 
chemical bonding. Fig.[TJb-c) illustrates this charge density 




FIG. 2: (color online) Ab initio COHP analysis of the orbitals participating to the bonding of the oxygen atom in epoxy position with 
its neighboring carbon atoms. Region of interest between — 1 eV and 1 eV with respect to the Fermi level. Positive and negative values 
indicate bonding and anti-bonding interactions, respectively. 



difference p defined as follows: 



where p° and pS m P h represent the charge densities of free- 
standing oxygen and graphene respectively. p tot is the charge 
density of the total system in its bound state. 

In Fig.|TJb), the charge accumulation in red and the charge 
depletion in blue indicate a covalent bonding between the 
carbon and the oxygen atoms. Oxygen is known to exhibit 
an acceptor behavior (attracts electrons on the s and p y or- 
bitals). Its p y orbital does not participate in the bonding 
as the electronic cloud around the oxygen keeps its conven- 
tional p-orbital form. 

In Fig. [TJc) the accumulation of electrons is rendered in 
red and green and the depletion in blue. The typical w- 
orbitals of sp 2 graphene are broken by the epoxy bonds. 
Indeed, this bonding causes a charge depletion (blue re- 
gion) in the ir electron cloud located above and under the 
d bond associated with the underlying carbon atoms. The 



latter a bond thus encounters a slight charge accumulation 
(green region). Additionally, the electronic charge trans- 
fer calculated with a HirshfelcP^l ( VoronoP^) integration 
adds up to —0.256 |e| (—0.266 |e|) of charge transfer to- 
wards the oxygen (red and green region close to the oxygen 
atom), which indicates that the bonding is also slightly ionic. 
These two integration techniques yield better results than 
the frequently used Mulliken integration since they are basis 
independent 77 . The oxygen bound in epoxy position is found 
to perturb graphene, but not in a purely sp 3 -hybridization 
way as in the hydrogenation case. 

A Crystal Orbital Hamilton Population (COHP) study^l221 
partitions the band structure energy in terms of orbital pair 
contributions. Positive density values represent the bond- 
ing states, while negative values describe anti-bonding states 
when plotting the conventional -COHP. Fig. [2] presents a 
COHP study on the orbitals participating to the bonding 
between the oxygen atom in epoxy position and its neigh- 
boring carbon atom. In the following analysis, conclusions 
are drawn for the region of interest for transport properties 
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FIG. 3: (color online) Electronic band structure of a 5 x 5 graphene 
supercell containing a single epoxy group. The DFT (red-solid) 
and TB (blue-dashed) band structures (a) are described along high- 
symmetry paths as described in the Brillouin zone (b). Nomencla- 
ture of extracted TB parameters is illustrated in (c). One possible 
orientation of the oxygen atom in epoxy position is described here- 
with. The two other inequivalent orientations are obtained by rotat- 
ing the oxygen atom in epoxy position by 120° around the central 
carbon atom in (d). 



[—1 eV; 1 eV] only; given bonds may have different bonding 
or antibonding behaviors away from the Fermi energy. 

One first notes that the contributions of the s (a) and p z 
(c) orbital of oxygen bonding with a neighboring carbon are 
analog. Both these orbitals have a dominant antibonding 
contribution with the p x and s orbital of carbon at the right of 
the Fermi energy. The COHP study between both first neigh- 
boring carbon atoms [Fig.[2](a), inset] shows that these two 
orbitals of carbon, also responsible for the a bond hybridiza- 
tion, strongly bind at the same energies (E ~ 0.5 eV and 
1.2 eV), with analog contributions, typical for the er-like hy- 
bridization between carbon atoms. These four orbitals (O p ,, 
O s , C Pt and C s ) thus form a hybridized electronic cloud 
with bonding contributions between the two carbon atoms 
and antibonding contributions between oxygen and carbon. 
This result is in agreement with the charge density rearrange- 
ment observed in Fig.[T] The electronic charge depletion of 
the 7T orbitals observed in Fig. [T] can be rationalized with a 
COHP study for a larger energy window (not shown here). 
The 7r electrons are partly drawn into the a bond between the 
carbon atoms. The remaining p z electrons of carbon [blue 




FIG. 4: (color online) DFT STM image of empty states, ie. Local 
DOS for one oxygen atom in epoxy postion, integrated between 
and 0.5 eV. Both sublattices are affected equally preserving the 
pseudospin symmetry in graphene. 



lines in (a), (b) and (c)] interact mainly with the p x , p z and 
s orbitals of oxygen. Finally, the p y orbital of oxygen does 
not interact with carbon in this energy window. 

Consequently, a TB model with first-neighbor interactions 
including both p x and hybridized slp z orbitals of oxygen 
with the p z and s/p x orbitals of carbon should be sufficient 
to accurately model the effect of epoxy groups on graphene. 
The combined contributions of the three orbitals of carbon 
binding with oxygen is renormalized to only one orbital in a 
7r-like model. The matrix elements of the TB Hamiltonian 
are given by: 

ij i 

At first, a \/3 x \fi R 30° supercell with one epoxy atom was 
simulated using DFT to extract a band structure with limited 
folding of the Brillouin zone (not shown here). The bands 
near the Fermi energy are nicely fitted using the following 
TB parameters [nomenclature, see Fig.(3jc)]: e x — —2.5 eV, 
e z = -1.0 eV, e 1 = 1.5 eV, j x = I.870, j z = -1.57 and 
71 = 0.0 eV with 70 = -2.6 eV. 

Finally, to confirm the validity of the model, these TB pa- 
rameters are used to generate the band structure of a 5 x 5 
supercell containing a single epoxy oxygen, which is super- 
imposed with its DFT counterpart [see Fig. [3ja)] . The elec- 
tronic path chosen to plot the band structure contains all in- 
equivalent high-symmetry segments in the 2D Brillouin zone 
of a 5 x 5 graphene supercell 80 [see Fig.[3jb)]. Note the band 
crossing at the Fermi energy is shifted away from the K2 
point. The localized fiat band (visible around the T point) 
appearing in the DFT band structure around -2.5 eV is miss- 
ing in the TB band structure. This originates from a strong 
interaction between the p y orbital of oxygen with the p y or- 
bital of carbon. Both these orbitals are missing in the TB 
model. These band structures were calculated for one ori- 
entation of the epoxy group on graphene[see Fig. [3jd)]. For 
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the two other possible orientations of the epoxy oxygen the 
crossings occur close to K\ or K3 for symmetry reasons. 
In conclusion, our TB model seems to be sufficient to ac- 
curately model random positions and random orientations 
of impurities, as long as these epoxy oxygens do not inter- 
act with each other which is assumed to be satisfied for the 
range of concentrations of impurities considered here. 

Finally, in Fig. |4] the twofold D 2 h symmetry in the sim- 
ulated STM image obtained by integrating the local DOS 
(LDOS) proves the analogy with a double impurity defect, 
by comparison with the LDOS in Fig. 2 of Ref. |60|. The 
LDOS of empty states is spatially integrated between and 
0.5 eV. A similar pattern (not shown here) is obtained for 
hole carriers, by integration between and and —0.5 eV. 



B. Kubo Formalism 

Transport properties for large mesoscopic-sized systems 
can be simulated efficiently using an order-N method based 
on the Kubo formalisrrF^H^I Assuming the electronic trans- 
port in the system is isotropic for the in plane x and y direc- 
tions, the 2D diffusion coefficient D(t) is obtained by 



D(t)=D x (t) + D y (t)=2D x (t) 



(3) 



Within this formalism, the diffusion coefficient D x (t) in the 
transport direction x is calculated at each time step using 



D x (t) = 



AX 2 {t) 



(4) 



where 



AX 2 (E,t) 



Tr 


5(E - H) 


X(t)- 


1(0) 


2" 


Tr 


S(E - H) 





(5) 



used for metals, which considers that the oscillation of the 
recursion coefficients is rapidly damped with the number 
of recursion steps. We checked that for the recursion step 
n = 500 the damping is sufficient, although a very small 
remnant oscillation caused by the small energy gap at high 
is observed. This is quantitatively correct at low energies 
and qualitatively sufficient at the border of the energy spec- 
trum where the energy gaps occur, by comparison with other 
more sophisticated termination methods. 

From the diffusion coefficient, the mean free path £ e (E) 
and the semi-classical conductivity <r sc (E) can be calculated 
using respectively: 



D mdx (E) 
2v(E) 



and 



a sc (E) = -e 2 p(E)D™*(E) 



(7) 



(8) 



where v(E) is the charge carrier velocity at energy E, D max 
the maximum value of D(t), e the electronic charge, and 
p(E) the DOS at energy E. The semi-classical Kubo- 
Greenwood conductivity er sc can be compared to the Drude 
approximation close to the Dirac point: 



4e 2 k(E)£ e (E) 



(9) 



h 2 

where E — fivpk with vp the Fermi velocity close to the 
Dirac point (d cc ps 1.42 A): 

37o d c 



v F 



2h 



1 x 10 D ms~ 



(10) 



The limitations of the Drude approximation have recently 
been discussed and put into context in a Review paper [90] 
and the limitations of the Born approximation in the Boltz- 
mann theory of conductivity have been analyzed in Ref. |9"T1 . 



IV. RESULTS 



where X{t) is the position operator in Heisenberg represen- 
tation at time t: 



X(t) = U\t)X{Q)U{t) 



(6) 



and where U{t) — e - %Ht / h is the time-evolution operator. 

The trace, which is a sum over wavepackets initially lo- 
calized on each orbital of the system, is replaced by an ini- 
tial state with a random phase on each orbital of the sys- 
tem. Taking the average of ten initial random phase states 
already yields very satisfactory results on the smoothness of 
the curves. This greatly reduces computation time. U(t) can 
be expanded using Chebyshev polynomials to allow for the 
mandatory order-N method to achieve reasonable computa- 
tion time for systems containing millions of orbitals. Both 
the numerator and denominator in Eq. |5]l are calculated us- 
ing the Lanczos recursion scheme thanks to continued frac- 
tions expansions. The termination term is the one usually 



To begin with, we discuss the effect of various oxygen 
concentrations on the density of states (DOS) in compari- 
son with pristine graphene. Then, the results obtained for 
the diffusion coefficient D(t) are analyzed. The other trans- 
port quantities are calculated within the Kubo formalism as 
introduced in the previous paragraph. Particular attention 
is given to the scaling behavior of the Kubo conductivity. 
All TB calculations were performed on systems containing 
2560000 carbon atoms, which corresponds approximately to 
systems of 300 nm by 200 nm. 



A. DOS and energy shift 



The evolution of the DOS, calculated from 



Tr 



S(E-H) 



with increasing impurity density is re- 



ported on Fig. [5] Although DFT calculations on similar 
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FIG. 5: (color online) Main frame: DOS for various impurity den- 
sities ranging from 0.05 to 4.42 %. The minimum of DOS shifts 
slightly with the addend concentration due to the changes in hop- 
ping parameters and on-site energies in the TB model. Side panels: 
(5-like peaks corresponding to impurity bands at lower (left) and 
higher (right) energies. 
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FIG. 6: (color online) Main frame: Realignment of the minimum of 
DOS and charge neutrality point (CNP: position obtained by DOS 
integration) at eV. Inset: linear increase of energy shift Ae n with 
increasing concentrations of epoxy groups (x in %). Rigid band 
theorem (see text) implies an impurity induced potential of ~ —3.7 



concentrations exhibit a shift of the Fermi energy compared 
to the pristine case (not shown here), no shifts are applied 
yet to point out the similarity with other studies^EI cov- 
ering the effect of ideal impurities (z'e. the influence of the 
change of respectively the on-site energy and the hopping 
parameters) on the DOS of pristine graphene. 

A shift of the minimum of the DOS is observed with in- 
creasing impurity concentrations. Even though formally the 
rigid band theorem cannot be used for atoms that do not 
have the same valence (ie. carbon and oxygenj^ll the shift 
is found to be linear with increasing concentrations {x) and 
the second order corrections 0(2) can thus be neglected [see 
Fig.[6j(inset)]: 



Ae„ = xU + 0(2) 



(11) 



where U is the local potential induced by the epoxy defect. 
A linear fit of Ae„ versus x implies a value of U equal to 
- -3.7 eV. 

In addition, the Van Hove Singularities (VHS) are 
smoothened out and decreased in amplitude with increasing 
concentrations of epoxi de groups, in agreement with previ- 
ous observation :! 82 ! 92 ! 93 -!. Also, a small increase of the DOS 
appears at the minimum of the DOS with increasing impu- 
rity concentration. 

The bumps in densities of states on the left and on the 
right of the minimum of DOS (Fig. [5] middle panel) cor- 
respond to the resonant energies between the oxygen atoms 
in epoxy position and the graphene sheet 60 , while the (5-like 
peaks in the side panels at high energies correspond to flat 
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impurity bands. As discussed in Section III A the localized 



FIG. 7: (color online) Projected densities of states. The p x (p z ) 
orbital of oxygen mainly contributes on the left (right) side of 
the Fermi energy. Negligible contributions are predicted for first- 
nearest neighboring carbon atoms. Second carbon nearest neigh- 
bors do contribute to the resonant energy bump with oxygen. 



state close to the left VHS caused by the strong interaction 
between both p y orbitals of oxygen and carbon is missing in 
this simplified TB model. 

Finally, oxygen in epoxy position triggers a shift of the 
Fermi energy which compensates the shift of the minimum 
of DOS discussed above (see Fig. 15} . This Fermi energy 
is obtained by integrating the TB DOS and counting the 
number of electrons present in the system. In the next Sec- 
tions, transport calculations will implicitly include this re- 
alignment of minimum of DOS with the Fermi energy, thus 
locating the charge neutrality point (CNP) at E — eV. 

In Fig. [7] a Projected Densities of States (PDOS) evidence 
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FIG. 8: (color online) Time evolution of diffusion coefficient (nor- 
malized) at E = 0.5 eV for different impurity concentrations. 



FIG. 9: (color online) Normalized diffusion coefficient versus time 
for x = 1.77 % at different energies. 



that the bump on the right side of the Fermi level originates 
from the p z orbital and the one on the left from the p x or- 
bital of oxygen, in agreement with the previous COHP DFT 
study. Fig.[7]also indicates that the first-nearest neighboring 
carbon atoms do not contribute to the total DOS in contrast 
to the second nearest neighboring carbon atoms. The PDOS 
suggest the oxygen atom atttacts most of the electtonic den- 
sity and thus weakens the density on the first-neighboring 
carbon. Such analysis agrees with existing literatu re 95 * 96 !. 



B. Diffusion coefficient and transport regimes 

The diffusion coefficient D (t) inherently contains all the 
information needed to calculate transport properties [see 
Eqs. ([7]) and Its time evolution or dynamics also clar- 
ifies the dominant transport regime at the considered time 
scale (ie, ballistic, diffusive or localized). Fig. [8] illustrates 
the typical behaviors of the normalized diffusion coefficients 



D nom (t) 



D(t) 

JJmax 



(12) 



for an energy E — 0.5 eV. As expected, the conduction in 
pure graphene is simply ballistic. The smallest simulated 
impurity concentration (ie. 0.01 %) does not reach its max- 
imum for the diffusion coefficient for the total elapsed time 
considered here (approx. 3800 fs). The slope of D(t) gives 
access to a v 2 (E) value still in close agreement with the the- 
oretical (analytical) value vp for energies E close to Ep. 
For instance, at eV and for 0.01 %, v(E) = 2.HA70//1, 
while vp = 2.13A7o/7i. For intermediate densities (0.05 
and 0.1 %) D n0lm (t) evolve from the ballistic regime to the 
diffusive regime for the simulated times. At higher concen- 
trations (0.49 %) D norm (t) reach a diffusive regime sharply, 



followed by a clear decrease with time, signature of quantum 
interferences leading the charge carrier localization. 

Fig. [9] presents diffusion coefficients D noim (t) of 1.77 % 
of impurities at different energies. Localization effects fol- 
low an asymmetric behavior between electrons and holes. 
The predominant resonant peak at the right of the charge 
neutrality point (CNP) causes much stronger localization ef- 
fects than the weaker and more smeared out peak at the left 
of the CNP (in analogy with the PDOS in Fig. [7). Such 
asymmetty is of crucial importance to synthesize possible 
electronic devices made of functionalized graphene requir- 
ing an efficient switch between a conducting and an insu- 
lating behavioi'22122]. Localization effects are more signifi- 
cant around 0.75 eV. One finally notes the peculiar behavior 
of the diffusion coefficient exactly at the CNP. Numerically, 
this point is more problematic to simulate since the density 
of charge carriers may be very scarce. Nevertheless, our re- 
sults suggest that the saturation limit of the Diffusion coef- 
ficient is reached for longer simulation time and then turns 
into a moderate decrease characteristic of quantum effects. 

This discussion on the diffusion coefficient points out 
two different regimes which are approached separately in 
the remaining Sections. Firstly, the semi-classical quanti- 
ties, which neglect quantum effects, are discussed in Section 



IV C Secondly, the localization regime is analyzed in Sec- 
tion llV PI 



C. Semi-classical regime 

1. Analysis of elastic mean free paths 

Using Eq. Q, the mean free path l e is calculated using 
the diffusion coefficient and plotted in Fig. 10 for energies 



between the two VHS at —2.6 eV and 2.6 eV. For the con- 
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FIG. 10: (color online) Energy dependent mean free paths for 
different impurity concentrations. Inset: mean free path versus 
hole (dashed/squares) and electron (solid/circles) carrier density for 
epoxy concentrations of 0.49 % and 0.95 %. 



FIG. 1 1 : (color online) Ratio of mean free paths at two selected 
impurity densities. £ e (xi) /£ e (x2) (solid lines) and xz/xi (dashed 
lines) with x 2 = 0.49 %. 



sidered elapsed times, the diffusive regime is not reached 
at every energy for the smallest impurity concentrations (ie. 
fii < 0.1 %, see Fig. [8] upper panel). Therefore, the mean 
free path can not be estimated for these small concentra- 
tions. For concentrations larger than 0.5 %, the diffusive 
regime is reached within the entire energy window. The dip 
in the mean free path at the right of the CNP (and in a lesser 
extent, at the left of the CNP), indicating larger scattering 
effects, shifts away from the CNP and becomes smoother 
with increasing concentration of oxygen atoms in epoxy po- 
sition. Such dips correspond to the resonance peaks found 
in the DOS which are induced by the oxygen. The inset 
of Fig. 10 confirms the predicted asymmetry, affecting the 
electrons more strongly. The largest concentrations scatter 
more uniformly across the entire energy window. In addi- 
tion, the evolution of the mean free path with impurity con- 
centration follows a simple scaling law as expected from a 
Fermi golden rule (Fig. [TT|: 

(■e{X\) _ 2-2 
e e (x 2 ) Xl 



2. Mobility 



(13) 



In Fig. 12 (mainframe) the mobility of the charge carriers 
is estimated theoretically using: 



H(E) 



n e 



(14) 



Scattering effects are affecting the electron mobility more 
strongly than the hole mobility. This asymmetry is reduced 
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for the largest impurity concentrations [see Fig. 12 (inset)] 



FIG. 12: (color online) Main frame: electron and hole mobility for 
usual experimental carrier densities. More severe scattering effects 
for negative charge carriers cause an asymmetry in mobility com- 
pared to the hole mobility. Inset: same mobilities with respect to 
the gate voltage V g . 



Experimentalists usually consider the absolute value of 
the mobility as a key quantity to characterize samples and 
corresponding inherent disorder. Temperature breaks the 
phase coherence of electrons along the scattering path and 
generally reduces quantum interference effects. Accord- 
ingly, the use of the semi-classical conductivity in evalua- 
tion of fi(E) (Eq. 14 1 is a reasonable approximation to an- 
alyze the experimental data. On the same basis, computed 
semi-classical conductivities are expected to be more valu- 
able for comparison with conductivities measured experi- 
mentally at room temperature. One may argue that in such 
a non-zero temperature environment, electron-phonon cou- 
pling may play also a significant role. However, inelastic 
scattering lengths due to electron-phonon coupling are ex- 
tremely long in graphene and may thus be disregarded too, 
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1000 




FIG. 13: (color online) Comparison between the numerical Kubo 
conductivity <t sc (solid lines) and the Drude approximation ctd 
(dashed lines) for different impurity concentrations. Minimum the- 
oretical value 2/ttGo is plotted in horizontal dashed green line. 



stronger classical scattering effects on the same energy win- 
dow (see Fig.|9|. This will be emphasized in Section IV D 
by studying the evolution of the conductivity with time deep 
into the diffusive regime. 



4. Short-range vs long-range scattering 

The nature of scattering range induced by oxygen atoms 
placed in epoxy position can be discussed using the scaling 
properties of semi-classical quantities. Both the scattering 
time and the conductivity are here briefly outlined^. 

The elastic scattering time is defined by 



t(E) 



te(E) 

v(E)- 



(15) 



Using the dispersion relation E = hvpk, t(E) can be es- 
timated in terms of the Fermi wave vector k. According to 
Nomura et al. |99| who also used the Kubo-Greenwood ap- 
proach, following conclusions apply. Long range (lr) and 
short range (sr) disorder scattering times scale respectively 
as follows: 



at least as a first approximation. 



3. Numerical Kubo conductivity and Drude approximation 

Fig. [13] compares the semi-classical value of the Kubo- 
Greenwood conductivity cr sc (solid lines) with the Drude 
conductivity cjd (dotted lines), extracted from Eq. |8]l and 
(|9]l respectively. The Drude approximation seems to be valid 
only for small impurity concentrations in the energy window 
[—1 eV; 1 eV]. For larger densities, the conductivities are 
underestimated by a factor of two for energies close to the 
CNR Moreover, cr D is ill-defined at the CNP [see Eq. |[9)]. 
Indeed, at the CNP, k = and co nsequently, ctd vanishes. 
However, theoretical reports extensively on a 

minimum value of the semi-classical conductivity equal to 
2/ttGq = 4e 2 /irh (when neglecting quantum interferences). 
To avoid the singularity at eV, the mean free path could be 
calculated analytically, including its k dependence explic- 
itly. As illustrated by the calculations for larger impurity 
concentrations, the variations of the DOS with disorder have 
to be included, as it is presently the case in the Kubo formula 
[Eq. ([HJ]. These variations account for single and multiple 
scattering events mentioned in Ref. 11031 . Up to some nu- 
merical discrepancy, the semi-classical minimum value of 
conductivity is confirmed. 

Additionally, an asymmetry exists between the electron 
and hole conductivities. This effect is weakened but not 
cancelled with increasing impurity concentrations. This is 
similar to the asymmetry already observed for both charge 
carriers in their respective mean free paths and mobilities. 
Albeit physically different, the stronger quantum localiza- 
tion effects on the electron side are directly linked to the 



Tir cx k and r sl - oc — 
k 



(16) 



Following such criterion, our data with corresponding nu- 
merical fits (Fig. [T4) i clearly indicate a short-range scatter- 
ing behavior of oxygen in epoxy position. Such short-range 
scattering time should diverge at the CNR This does not hap- 
pen within the Kubo formalism since the DOS remains finite 
close to the CNP, in contrast with the prediction of the Drude 
approximation. 

A similar analysis based on the scaling of conductivity 
is not straightforward. According to Ref. |99l , long range 
conductivity should scale linearly with n, while short range 
conductivity should display a non-linear behavior, approach- 
ing the constant Boltzmann (or Drude) conductivity g-q for 
\Ep\ ^> Ti/t. In Fig. 13 the behavior close to the Dirac 
point behaves differently depending on the impurity concen- 
tration. For smaller impurity concentrations, the conduc- 
tivity tends to decrease for energies corresponding to rea- 
sonable carrier concentrations (up to 10 13 cm~ 2 ), while it 
increases slightly for larger impurity concentrations. The 
Drude conductivity never reaches a constant plateau for the 
whole of the energy E or carrier n range. The short-comings 
of this conductivity have already been pointed out. 

We note that this subject is still debated, as other 
model d 41 ! 91 ! 100 ! 101 ! predict a dominant linear dependency 
(with logarithmic corrections) for both long range and short 
range conductivity. 

Additionally, it has been observed experimentally^! that 
a quasi-ballistic regime exhibits a non-linear behavior with 
n, while a more disordered system scales linearly with n. 
Comparison with experiment becomes particularly tricky as 
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FIG. 14: (color online). Scattering times in function of Fermi wave 
vector k for different oxygen densities. Numerical fits for each data 
set based on a regression algorithm with fitting parameters A and 
B (black) to account for the generalized inverse function. Hole and 
electron scattering times were fitted separately. 



most of previously mentioned analytical models disregard 
important multiple scattering effects on the computed con- 
ductivity. 

Another important remark is that most theoretical predic- 
tions have been derived assuming restrictions on disorder 
models which are partly invalidated in the present study. 
Indeed, the epoxy defects have been derived from accu- 
rate first-principle calculations, and the resulting TB model 
brings more realism and generality when compared to sim- 
plified academic models. As a matter of fact, the DOS 
(Fig. [6]l evidences resonant energy bumps, driven by ran- 
domly distributed oxygen, which could cause squared loga- 
rithmic corrections^. Our data cannot really be accurately 
fitted to obtain the different corrections to the scaling in this 
context. 

This dominant short-range disorder is at the root of the 
quantum effects presented in the remaining Sections. 



D. Evolution of the Kubo conductivity with time scale (or 
length) 

The semi-classical expression of the Kubo-Greenwood 
conductivity [Eq. (7)] restricts the transport to the diffusive 
regime, ie. when suppressing quantum interferences. To fol- 
low the time evolution spreading of the quantum wave pack- 
ets, the expression of D x (t) in Eq. (3) should be replaced as 
follows: 



D x (t) = 



d(AX 2 {t)) 
dt 



FIG. 15: (color online) Kubo conductivities at different elapsed 
times of wave packet evolution for 4.42 % of impurities. 



When replacing £) max by the new expression of D(t) [Eq. 
Q] in Eq. (jHJ, the Kubo conductivities are obtained at dif- 
ferent time scales as depicted in Fig. [15] for a 4.42 % im- 
purity density. The time-evolution of the conductivity is 
found not to be uniform over the energy spectrum, indicat- 
ing there might be different transport regimes depending on 
the charge carrier energies and impurity concentrations for a 
given length scale. 

According to the scaling theory of localization 36 , there are 
two possible behaviors for the conductivity corresponding 
to the weak and strong localization regimes which read as 
follows: 



a{L) = a\ 



( Lit) 
hir 2 V y/2£ e 



In 



cr(L) ~ exp 



L(t) 



(18) 



(19) 



where the localization length £ gives an estimation of the dis- 
tance covered by a charge carrier before it is totally trapped 
due to this multiple scattering effects. The diffusion length 
is defined as 



L(t) = 2 y / 2AX 2 (t). 



(20) 



This definition of L is reasonable when saturation of the 
diffusion coefficient has been reached. The extra fac- 
tor y/2 in Eq. 



18 



compared to the correction obtained by 
Lee et al. [36| comes from a different definition of D (t). 
Both numerical estimation of a(L) (symbols) and analyti- 

2 /hn 2 In fromEq. ( 18 1 (solid lines) 

The numerical part con- 
tains small jiggling caused by the very sensitive derivation 



cal 



are plotted in Fig.[T6]for L > L 



(17) inEq. ( 17 i 
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FIG. 16: (color online) Kubo conductivities for L > L max , which 
include weak localization corrections to the semi-classical conduc- 
tivity, for different impurity densities at energy 0.8 eV. The numer- 
ically estimated conductivity u(L) (obtained using Eqs 



17 



and 



20 



symbols) contains numerical jiggling caused by the very sensitive 
derivation in Eq. ^] The conductivity obtained using Eq. j 1 8fr is 
plotted in solid lines. Only one out of five points are plotted in the 
inset for clarity reasons. 



The corrections to the semi-classical conductivity in the 
low impurity limit (mainframe) follows the logarithmic be- 
havior, from which an estimation of £ can be deduced. £ cor- 
responds to the length where the cooperon corrections equal 
the semi-classical conductivity 36 . Starting from Eq. 
localization length is thus estimated with 



18 1, the 



exp 



Go 



(21) 



with the computed i e and cr sc and corres ponds to the 2D gen- 
eralization of the Thouless relationship 104 ! 105 !. These local- 
ization lengths are plotted in Fig.flTpSl 

For larger concentrations, the cooperon corrections to the 
semi-classical conductivity seem to saturate and depart from 



the perfect logarithmic behavior (Fig. 16 inset). The cor- 
rections obtained numerically become smaller than what is 
predicted due to a transition to the strongly localized regime 
following an evanescent exponential behavior. This can be 
rationalized invoking the Ioffe and Regel criterion 107 which 
states that the appearance of strong localization becomes sig- 
nificant for impurity concentrations satisfying kp£ e = 1. 
Such a criterion implies a mean free path of the order of the 
interatomic distance, which is actually the case for 4.42 % of 
impurities where the mean free path is estimated to be below 
3 A for the energy window from 0.5 eV to 1.0 eV. 

In Fig. [18] by fitting the exponential behavior of Eq. ( fT9) , 
values for £ equal to 8.4 and 4.8 nm are estimated for 3.22 
% and 4.42 % of impurities respectively at an energy of 0.8 
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FIG. 17: (color online) Localization lengths estimated using 
Eq. {XT} for different impurity concentrations. 
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FIG. 18: Exponential fits [Eq. \19) ] to estimate localization lengths 
£ in the strongly localized regime at energy 0.8 eV. A first fit 
(dashed lines) for the whole of the available data allows us to esti- 
mate a length L for which wavepackets are localized. A second fit 
(solid lines) for values larger than L gives us new estimates for £. 



eV (dashed lines). Refitting cr{L) for the region at the right 
of these values (solid lines), we obtain convincing exponen- 
tial decays and more accurate estimates for £ equal to 11.2 
and 5.3 nm, respectively. Both these estimates and the ones 
obtained by Eq. pT| ) in Fig.|17|are of the same order of mag- 
nitude, thus validating our results. Experimentalists however 
often use the Drude approximation ajj instead of the correct 
semi-classical conductivity er sc in Eq. (21 1, The inaccuracy 



of the Drude approximation for largest impurity concentra- 
tions causes the localization length to be underestimated by 
an order of magnitude. 
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V. CONCLUSIONS 

In this paper, the quantum transport properties of chem- 
ically damaged two-dimensional graphene based structures 
have been investigated. Using the Kubo-Greenwood trans- 
port framework, and by means of an efficient order N numer- 
ical implementation, mesoscopic transport features in disor- 
dered graphene have been explored in details, with impuri- 
ties (adsorbed oxygen-driven epoxide defects) described by 
local tight-binding parameters, deduced from first-principles 
calculations. 

In addition to the numerical calculation of the energy- 
dependent elastic mean free path driven by a given epox- 
ide density, quantum localization effects have been analyzed 
from the weak to the strong (Anderson) localization regimes. 
By applying the conventional scaling theory of localization, 
the 2D-localization lengths have been evaluated from the 
scaling behavior of the Kubo conductivity, and contrasted to 
the prediction deduced from the cooperon correction to the 
conductivity (which relates £ to the elastic mean free path 
and semi-classical conductivity). A very reasonable agree- 
ment has been obtained, pinpointing further towards a strong 
energy-dependence of all transport length scales. 

By combining the ab-initio approach for the description 



of the defects structure and local energetics with an efficient 
and exact quantum transport methodology implemented on 
tight-binding models, our general theoretical framework 
provides a solid foundation and tool to understand the origin 
of complex transport phenomena in strongly disordered and 
chemically complex graphene-based nanostructures. The 
extension of our study to any other kinds of defects (topo- 
logical, chemical, etc) and other types of two-dimensional 
structures is straightforward. 
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